class MDS:

    def euclidean_dist(self, x1, x2):
        return np.linalg.norm(x1 - x2)

    def get_dist_matrix(self, data):
        """
        num = data.shape[0]
        D = np.zeros((num,num),dtype=float)
        for i in range(num):
            for j in range(num):
                dist = self.euclidean_dist(data[i], data[j])
                D[i,j] = dist
        return D
        """
        expand_ = data[:, np.newaxis, :]
        repeat1 = np.repeat(expand_, data.shape[0], axis=1)
        repeat2 = np.swapaxes(repeat1, 0, 1)
        D = np.linalg.norm(repeat1 - repeat2, ord=2, axis=-1, keepdims=True).squeeze(-1)
        return D

    def get_maxtrix_B(self, D):
        DD = np.square(D)
        sum_ = np.sum(DD, axis=1) / D.shape[0]
        Di = np.repeat(sum_[:, np.newaxis], D.shape[0], axis=1)
        Dj = np.repeat(sum_[np.newaxis, :], D.shape[0], axis=0)
        Dij = np.sum(DD) / ((D.shape[0]) ** 2) * np.ones([D.shape[0], D.shape[0]])
        B = (Di + Dj - DD - Dij) / 2
        return B

    def fit(self, data, low_dim=2):
        D = self.get_dist_matrix(data)
        B = self.get_maxtrix_B(D)
        eig_values, eig_vects = np.linalg.eigh(B)
        eig_values_sort = np.argsort(-eig_values)
        values_sort = eig_values[eig_values_sort]
        vects_sort = eig_vects[:, eig_values_sort]
        top_value_diag = np.diag(values_sort[0:low_dim])
        top_vects = vects_sort[:, 0:low_dim]
        Z = np.dot(np.sqrt(top_value_diag), top_vects.T).T
        return Z

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_swiss_roll

X, t = make_swiss_roll(n_samples=5000, noise=0.2, random_state=42) # X为坐标 t为颜色

from MDS import MDS

model = MDS()
Z = model.fit(X,low_dim=2)
axes = [-11.5, 14, -2, 23, -12, 15]
fig = plt.figure()
ax = fig.add_subplot(121, projection='3d')
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=t, cmap=plt.cm.hot)
ax.view_init(10, 60)
ax.set_xlabel("$x$", fontsize=18)
ax.set_ylabel("$y$", fontsize=18)
ax.set_zlabel("$z$", fontsize=18)
ax.set_xlim(axes[0:2])
ax.set_ylim(axes[2:4])
ax.set_zlim(axes[4:6])
plt.title('3D swiss roll')
ax2 = fig.add_subplot(122)
ax2.scatter(Z[:, 0], Z[:, 1], c=t, cmap=plt.cm.hot)
ax2.set_xlabel("$x$", fontsize=18)
ax2.set_ylabel("$y$", fontsize=18)
plt.title('after MDS')
plt.savefig("figure1.png")